/*
 * The copyright of these materials are subject to the
 * original authors.
 * 
 * This code is distributed freely.
 */
package org.nebulaframework.benchmark.scimark2;

/**
 * LU matrix factorization. (Based on TNT implementation.) Decomposes a matrix A
 * into a triangular lower triangular factor (L) and an upper triangular factor
 * (U) such that A = L*U. By convnetion, the main diagonal of L consists of 1's
 * so that L and U can be stored compactly in a NxN matrix.
 * 
 * 
 */
public class LU {
	/**
	 * Returns a <em>copy</em> of the compact LU factorization. (useful mainly
	 * for debugging.)
	 * 
	 * @return the compact LU factorization. The U factor is stored in the upper
	 *         triangular portion, and the L factor is stored in the lower
	 *         triangular portion. The main diagonal of L consists (by
	 *         convention) of ones, and is not explicitly stored.
	 */

	public static final double num_flops(int N) {
		// rougly 2/3*N^3

		double Nd = (double) N;

		return (2.0 * Nd * Nd * Nd / 3.0);
	}

	protected static double[] new_copy(double x[]) {
		int N = x.length;
		double T[] = new double[N];
		for (int i = 0; i < N; i++)
			T[i] = x[i];
		return T;
	}

	protected static double[][] new_copy(double A[][]) {
		int M = A.length;
		int N = A[0].length;

		double T[][] = new double[M][N];

		for (int i = 0; i < M; i++) {
			double Ti[] = T[i];
			double Ai[] = A[i];
			for (int j = 0; j < N; j++)
				Ti[j] = Ai[j];
		}

		return T;
	}

	public static int[] new_copy(int x[]) {
		int N = x.length;
		int T[] = new int[N];
		for (int i = 0; i < N; i++)
			T[i] = x[i];
		return T;
	}

	protected static final void insert_copy(double B[][], double A[][]) {
		int M = A.length;
		int N = A[0].length;

		int remainder = N & 3; // N mod 4;

		for (int i = 0; i < M; i++) {
			double Bi[] = B[i];
			double Ai[] = A[i];
			for (int j = 0; j < remainder; j++)
				Bi[j] = Ai[j];
			for (int j = remainder; j < N; j += 4) {
				Bi[j] = Ai[j];
				Bi[j + 1] = Ai[j + 1];
				Bi[j + 2] = Ai[j + 2];
				Bi[j + 3] = Ai[j + 3];
			}
		}

	}

	public double[][] getLU() {
		return new_copy(LU_);
	}

	/**
	 * Returns a <em>copy</em> of the pivot vector.
	 * 
	 * @return the pivot vector used in obtaining the LU factorzation.
	 *         Subsequent solutions must permute the right-hand side by this
	 *         vector.
	 * 
	 */
	public int[] getPivot() {
		return new_copy(pivot_);
	}

	/**
	 * Initalize LU factorization from matrix.
	 * 
	 * @param A
	 *            (in) the matrix to associate with this factorization.
	 */
	public LU(double A[][]) {
		int M = A.length;
		int N = A[0].length;

		// if ( LU_ == null || LU_.length != M || LU_[0].length != N)
		LU_ = new double[M][N];

		insert_copy(LU_, A);

		// if (pivot_.length != M)
		pivot_ = new int[M];

		factor(LU_, pivot_);
	}

	/**
	 * Solve a linear system, with pre-computed factorization.
	 * 
	 * @param b
	 *            (in) the right-hand side.
	 * @return solution vector.
	 */
	public double[] solve(double b[]) {
		double x[] = new_copy(b);

		solve(LU_, pivot_, x);
		return x;
	}

	/**
	 * LU factorization (in place).
	 * 
	 * @param A
	 *            (in/out) On input, the matrix to be factored. On output, the
	 *            compact LU factorization.
	 * 
	 * @param pivit
	 *            (out) The pivot vector records the reordering of the rows of A
	 *            during factorization.
	 * 
	 * @return 0, if OK, nozero value, othewise.
	 */
	public static int factor(double A[][], int pivot[]) {

		int N = A.length;
		int M = A[0].length;

		int minMN = Math.min(M, N);

		for (int j = 0; j < minMN; j++) {
			// find pivot in column j and test for singularity.

			int jp = j;

			double t = Math.abs(A[j][j]);
			for (int i = j + 1; i < M; i++) {
				double ab = Math.abs(A[i][j]);
				if (ab > t) {
					jp = i;
					t = ab;
				}
			}

			pivot[j] = jp;

			// jp now has the index of maximum element
			// of column j, below the diagonal

			if (A[jp][j] == 0)
				return 1; // factorization failed because of zero pivot

			if (jp != j) {
				// swap rows j and jp
				double tA[] = A[j];
				A[j] = A[jp];
				A[jp] = tA;
			}

			if (j < M - 1) // compute elements j+1:M of jth column
			{
				// note A(j,j), was A(jp,p) previously which was
				// guarranteed not to be zero (Label #1)
				//
				double recp = 1.0 / A[j][j];

				for (int k = j + 1; k < M; k++)
					A[k][j] *= recp;
			}

			if (j < minMN - 1) {
				// rank-1 update to trailing submatrix: E = E - x*y;
				//
				// E is the region A(j+1:M, j+1:N)
				// x is the column vector A(j+1:M,j)
				// y is row vector A(j,j+1:N)

				for (int ii = j + 1; ii < M; ii++) {
					double Aii[] = A[ii];
					double Aj[] = A[j];
					double AiiJ = Aii[j];
					for (int jj = j + 1; jj < N; jj++)
						Aii[jj] -= AiiJ * Aj[jj];

				}
			}
		}

		return 0;
	}

	/**
	 * Solve a linear system, using a prefactored matrix in LU form.
	 * 
	 * 
	 * @param LU
	 *            (in) the factored matrix in LU form.
	 * @param pivot
	 *            (in) the pivot vector which lists the reordering used during
	 *            the factorization stage.
	 * @param b
	 *            (in/out) On input, the right-hand side. On output, the
	 *            solution vector.
	 */
	public static void solve(double LU[][], int pvt[], double b[]) {
		int M = LU.length;
		int N = LU[0].length;
		int ii = 0;

		for (int i = 0; i < M; i++) {
			int ip = pvt[i];
			double sum = b[ip];

			b[ip] = b[i];
			if (ii == 0)
				for (int j = ii; j < i; j++)
					sum -= LU[i][j] * b[j];
			else if (sum == 0.0)
				ii = i;
			b[i] = sum;
		}

		for (int i = N - 1; i >= 0; i--) {
			double sum = b[i];
			for (int j = i + 1; j < N; j++)
				sum -= LU[i][j] * b[j];
			b[i] = sum / LU[i][i];
		}
	}

	private double LU_[][];
	private int pivot_[];
}
